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ACCURATE  EFFICIENT  EVALUATION  OF  BESSEL 


TRANSFORM;  PROGRAMS  AND  ERROR  ANALYSIS 


INTRODUCTION 


The  method  of  Filon  integration  for  Fourier  transforms  [1],  [2;  pages 
408-409],  [3;  pages  67-75],  [4;  page  400],  [5;  page  890],  [6;  pages  62-66], 


dx  exp(ia>x)  g ( x ) 


is  well  established  and  very  useful  for  accurate  numerical  work.  Instead  of 
the  standard  Simpson's  rule,  which  would  approximate  the  complete  integrand 
exp(i(jx)  g(x)  by  parabolas  over  abutting  pairs  of  panels,  Filon's  method 
approximates  only  the  function  g(x)  by  parabolas,  and  carries  out  the 
corresponding  integrals  in  (1)  analytical ly.  These  closed  form  integrals 
are  then  evaluated  with  computer  aid.  Since  the  exponential  in  (1)  is  being 
handled  exactly  for  all  oj,  the  hope  is  that  the  error  of  approximating  (1) 
by  means  of  Filon's  method  will  be  substantially  the  same,  for  larger  o>  as 
for  small  u  (where  all  the  error  arises  from  approximating  g(x)).  That  is, 
Filon's  method  is  expected  to  be  an  error  maintenance  procedure,  whereby  the 
absolute  error  does  not  increase  significantly  with  o>.  Certainly  that  is 
not  the  case  for  the  Trapezoidal  and  Simpson  rules,  where  significant 
aliasing  severely  limits  the  accuracy  of  the  results  for  larger  a>. 
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An  alternative  simpler  procedure  to  Filon's  method  for  Fourier 
transforms  is  to  approximate  g(x)  by  straight  lines  over  abutting  panels, 
and  again  to  evaluate  the  resultant  integrals  in  (1)  analytically  in  closed 
form.  This  (less-accurate)  procedure  is  documented  in  [8;  pages  418-419], 
for  example. 


Here,  we  will  extend  these  two  procedures  to  a  Bessel  transform  of  the 


form 


G(oj)  = 


^  dx  JQ(wx)  g(x)  , 
0 


(2) 


where  g(x)  is  an  arbitary  given  function,  and  Jq  is  the  zeroth-order 
Bessel  function.  One  of  the  major  differences  we  encounter,  relative  to 
filon's  method,  is  that  the  resultant  integrals  cannot  all  be  evaluated  in 
closed  form.  In  order  to  circumvent  this  problem,  we  use  a  combination  of  a 
downward  recursion  and  an  asymptotic  expansion,  which  are  limited  in 
accuracy  only  by  the  inherent  round-off  error  of  the  computer  utilized, 
thereby  obtaining  an  efficient  useful  procedure  for  numerical  evaluation  of 
the  pertinent  integrals  and  functions. 


To  give  a  physical  application  where  the  Bessel  transform  arises, 
consider  that  we  are  interested  in  two-dimensional  Fourier  transform 

-eo6 


ff 


dx  dy  exp(iux  +-  ivy)  f^fx.y) 


-06 

where  function  f  has  isotropic  behavior.  That  is,  suppose  the  dependence 
of  f  is  solely  on  the  distance  from  the  origin  of  coordinates: 


(3) 


f2( x.y) 


fi< 


(4) 


Then  (3)  becomes 


+*6 


^  dx  dy  exp(iux  +  ivy)  f^(|/x^  e  y2) 


00 


=  2ir  |  dr  Jo(cor)  r  f^r)  , 


(5) 


where  we  changed  to  cylindrical  coordinates  and  have  defined 


,  2  21/2 
CJ  =  (u  +  V  ) 


(6) 


Thus,  (5)  is  of  the  form  of  (2),  upon  identification  of  g(x)  as  x  f^x). 


Suppose  in  (3)  that  the  dependence  on  x,y  is  more  general  than  (4), 
namely  of  the  form 


which  allows  for  a  general  center  point  of  symmetry  x0.yQ.  as  well  as  a 
tilted  elliptical  shape.  Then  substitution  in  (3)  yields,  after  a 
cylindrical  coordinate  change,  the  result 


(1 


2wab 


exp(iuxQ  e  ivyQ) 


od 

|  dr  J(o)r)  r  f^r)  , 


(8) 


where  now 
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Again,  the  fundamental  Bessel  transform  of  the  form  of  (2)  results,  where 
g(x)  is  x  f^x). 


On  the  other  hand,  if  G(u)  is  specified  in  (2)  for  w  >  0,  the 
corresponding  solution  to  this  integral  equation  is 


f 


g(x)  =  x  \  dw  JQ(X(j)  o>  G(w)  , 


(10) 


which  is  again  a  Bessel  transform  of  the  form  of  (2). 


Thus,  we  have  presented  several  instances  where  the  transform  given  by 
(2)  is  of  interest  and  must  be  accomplished  accurately  for  large  as  well  as 
small  arguments  of  the  transform  variable  u. 


4 


LINEAR  APPROXIMATION 


The  integral  of  interest  here  is 

x 

r 

G(w)  =  ^  dx  Jq(u)X)  g(x)  ,  (11) 

where  left-end  point  x^  could  be  zero,  and  right-end  point  x^,  could  be 

taken  so  large  that  g(x)  is  essentially  zero  for  x  >  x  .  (If  x.  is 

negative,  the  values  of  g  could  be  folded  over  to  the  positive  x-axis,  using 

g(x)  +•  g(-x)  as  the  new  integrand,  since  Jq(wx)  is  even  in  x.)  We  break 

interval  x.  ,x  into  a  number  of  abutting  panels,  each  of  the  same  width  h, 

^  r* 

and  fit  g(x)  by  straight  lines  over  each  of  those  panels.  The  fits  for  the 
’eft-end  point  and  an  abutting  (internal)  point  xp  are  depicted  in  figure 
1,  where  it  is  temporarily  presumed  that  the  adjacent  sample  values  of 


Figure  1.  Linear  Approximat ions  to  g ( x ) 
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function  g(x)  are  zero;  this  allows  us  to  isolate  the  contribution  of  each 
sample  of  g(x)  to  tl^  total  desired  in  (11).  The  straight  lines  pass 
through  the  function  value  gn  =  g( x^)  at  sample  value  xn,  and  are  zero 
at  the  adjacent  sample  points,  h  is  the  sampling  increment  in  x  applied  to 
g(x).  The  situation  at  the  right  end  is  the  mirror  image  of  that  at  the 
left  end,  depicted  in  figure  1. 


If  u>  is  zero  in  (11),  the  approximation  afforded  to  the  integral  by 
means  of  figure  1  is  obviously 

S(0)  *  $  ,  *  Vl  *•••*  Vl  2  9rJ 


91  w 


for  u  0 


'  n 


which  is  just  the  Trapezoidal  rule.  For  u>  >  0,  considerably  more  effort  is 
required;  there  is  no  need  to  consider  u  <  0,  since  JQ(u>x)  is  even  in  u> 
Before  we  get  into  that  derivation,  we  must  introduce  some  auxiliary 
functions . 


SPECIAL  FUNCTION  DEFINITIONS 


Define  the  integral 


A(u) 


I 


dt  JQ(t) 


(13) 


6 


*  .*  **• 
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!v 


This  function  cannot  be  evaluated  in  closed  form;  a  table  of  A(u)  is 
available  in  [5;  pages  492-493].  On  the  other  hand,  the  integral 


J*  dt  t  JQ(t)  =  u  J1 (u) 
0 


3 

& 


is  immediately  available  by  use  of  [5;  9.1.30].  And  two  integrations  by 
parts,  coupled  with  (13),  yields  the  result 
u 

J  dt  t2  JQ(t)  =  u2  J^u)  +  u  JQ(u)  -  A(u)  . 


% 

•a 


We  will  also  find  use  for  the  auxiliary  functions 


B0(u)s  A(u)  -  u  JQ(u)  =  f  dt  t  (u  -  t)  JQ(t)  , 


u 

B^u)*  A(u)  -  J^u)  =  |  dt  (l  -  JQ(t) 


All  of  these  functions.  A,  Bq,  ,  are  zero  at  the  origin  and  are  odd. 
Numerical  evaluation  of  these  functions  is  considered  in  appendix  A. 
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ABUTTING  POINT 


For  an  abutting  (Internal)  point  in  the  interval  (x^.x^),  as 
depicted  on  the  right-hand  side  of  figure  1,  the  contribution  to  integral 
(11),  due  to  this  single  sample  point  gn  =  9(xp)»  is 


I 

n 


X  A 

J 

vh 


dx 


J  (coX) 
o' 


(i  +  y)  + 


x  +h 

J 


J  (wx) 
0 


(i  -  y)  , 


(18) 


where  we  have  defined 


We  now  assume  that  the  n-th  sample  point  xn  is  taken  such  that 


(19) 


xp  =  n  h  for  JL  <  n  <  r  .  (20) 

This  makes 


=  j  =  rational  .  ( 21 ) 

This  constitutes  a  restriction  on  ratio  xr/Xj  in  (11);  it  has  been  adopted 

here  in  order  to  minimize  the  number  of  calculations  of  the  Bessel  function 

J  later,  when  we  consider  the  multiple  values  of  u>  desired  for  (11). 
o 

(The  procedure  presented  here  can  be  extended  to  the  general  case  where  ^ 
is  arbitrary  and  xp  =  +  nh,  if  desired.)  If  x^  is  zero,  then  the 

choice  in  (20)  is  no  restriction  at  all. 
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APPROXIMATION  TO  INTEGRAL 


An  important  parameter  in  this  numerical  integration  procedure  is  the 


quantity 


8  =  tdh 


which  is  the  product  of  "radian  frequency"  u>  and  the  sampling  increment  h. 

As  we  shall  see,  values  of  0  near  ir  and  2n  will  constitute  points  of 
considerable  aliasing;  see  [4;  page  400]  for  a  discussion  of  the  Fourier 
transform  case. 


When  the  procedure  in  (18)— (19)  is  extended  to  include  the  left-end  and 
right-end  points  of  integral  (11),  and  the  various  integrals  evaluated  with 
the  help  of  (13)-(17),  the  total  approximation  is  given  by  appendix  B  in 
several  alternative  forms,  one  of  which  is  (B-7): 


<*>  G(«)  =  [igJl+1  -  a  *  D^]  B^e)  -  g^  J^e) 


[r  9r_1  -  (r  -  l)gj  B^re)  +  gr  J]  ( re)  + 


n[Vl  '  *  Vl]  Vne> 


n-i+1 


where 


9n  =  9(*n)  =  g(nh)  .  (24) 

Reasons  for  this  grouping  of  terms,  including  speed  of  execution  and  storage 
requ i rements ,  are  discussed  below. 


•„fcf  V.*'  '  •  • 
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SAMPLING  INCREMENT  FOR  w 

When  output  variable  <o  in  integral  (11)  is  restricted  to  multiples  of  a 
sampling  increment  A,  according  to 

=  kA  for  1  <  K-|  <  k  <  K2  ,  (25) 

then  ne  =  nkAh,  meaning  that  the  arguments  of  the  B^u)  function  in  (23) 
are  limited  to  integer  multiples  of  hA,  the  product  of  the  sampling  increment 
in  input  variable  x  and  the  sampling  increment  in  output  (transform) 
variable  ok  The  explicit  relationship  for  G(u>)  =  G(kA)  is  given  by 
specializing  (23)  to  the  values  (24),  thereby  obtaining 

i 

I 

kA  G(kA)  ~  (X  *■  ))9gJ  B]  CfkAh)  -  J^tikAh)  ► 

+  [r  gr-l  "  (r  “  1)9rl  (  rk Ah )  +  9r  ( rkAh)  * 

r-1 

+  ^  n[gnt-1  '  2gn  *  gn-l]  Bi(nkah)  for  k  --  1  •  (26) 

n=K+l 

i  COMPUTATION  TIME  CONSIDERATIONS 

I 

Thus,  we  need  evaluate  B^(u)  only  at  u  =  mAh,  where  m  is  an  integer. 
Furthermore,  not  all  values  of  integer  m  will  be  encountered  as  n  and  k  sweep 
out  their  respective  values  given  by  (20)  and  (25).  And  since  B ^ ( u ) , 
defined  in  (17)  and  (13),  is  the  most  time-consuming  aspect  of  the 
computation  of  (26),  it  behooves  us  not  to  compute  B^mAh)  at  values  of  m 
that  will  not  be  encountered,  and  not  to  recompute  B^mAh)  at  values  of  m 


W. 
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that  are  encountered  more  than  once.  This  latter  situation  arises  when  m  is 
highly  composite;  for  example,  m  =  12  =  4*3  =  6*2  =  12*1  could  be 
encountered  several  times  as  n  and  k  vary  in  (26). 

In  order  to  incorporate  this  time-saving  feature  into  the  Bessel 
integral  evaluations  required  by  (26),  the  values  of  B  (nkAh)  are  computed 
only  once  and  stored  in  a  one-dimensional  array  at  linear  location  m  -  nk. 
Unf ortunately,  this  speed-up  feature  is  achieved  at  the  expense  of 
considerable  storage,  for  if  n  and  k  range  up  to  N  and  K,  respectively,  the 
one-dimensional  storage  array  must  have  NK  cells,  of  which  most  are  empty 
when  N  and  K  are  large. 

when  N  and  K  are  so  large  that  storage  is  not  feasible,  such  as  when 
xr  in  (11)  is  large,  and  large  u  is  desired  in  (25),  then  the  alternative 
procedure  of  direct  brute-force  evaluation  of  (26)  for  B^(nkAh),  repeated 
as  often  as  necessary,  but  without  storage,  is  employed.  Recomputation  of 
B^(mAh)  for  some  m  values  occurs,  but  evaluation  at  unused  m  values  never 
does . 


Thus  we  have  two  alternatives  and  two  corresponding  programs  for  (26): 
one  faster  routine  which  may  require  considerable  storage,  and  a  slower 
procedure  utilizing  very  little  storage.  The  former  is  recommended  when 
feasible,  while  the  latter  furnishes  a  back-up  position.  Programs  for  both 
procedures  are  listed  in  appendix  B. 
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BEHAVIOR  FOR  SMALL  0 

When  ©  is  small,  differences  of  functions  w  th  similar  values  are 

required  in  (23),  as  may  be  observed  by  the  linear  u  dependence  on  the 

left-side.  The  appropriate  series  development  for  this  linear  approximation 

2 

approach  to  (11)  is  given  in  ( B-l 1 ) -( B-l 2) ,  through  order  0  .  Additional 
4  6 

terms  to  order  e  ,  e  can  be  derived  by  extending  the  approach  given 
there;  however,  an  easier  technique  will  be  developed  in  the  next  section. 


1  2 
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contribution  of  each  sample  value  gn  =  g(xn)  is  isolated,  by  temporarily 
presuming  that  the  neighboring  sample  values  of  g(x)  are  zero.  The  variable 
y  in  figure  2  is  again  the  normalized  quantity 


where  h  is  the  sampling  increment  in  x  applied  to  g(x). 

If  co  is  zero  in  (27),  the  approximation  afforded  to  the  integral  by 
means  of  figure  2  is 

G(0)  >5^  -  49<,1  *  2V?  ?9r-2  *  4Vl  *  9r]  = 

(29) 

=  3  [fCXf  )  f  4g(  x^i )  ?9(x#+2V- •  29(xr-2)  f  4g(xr-l}  f  g(xr0  ’ 

which  is  Simpson's  rule. 

APPROXIMATION  TO  INTEGRAL 

Since  Jq  in  (27)  is  even  in  <o,  we  only  need  to  consider  co  >  0  in  the 
following.  The  derivation  of  the  approximation  to  integral  (27),  by  means 
of  the  parabolic  fits  in  figure  2,  is  carried  out  in  appendix  C,  culminating 
in  ( C -1 0) -( C -1 2 )  : 


14 


Cm.X'JX-MjCd.TjlKnJ'tL 


lilMVAWJWWimV'A'A'JV'. 


;; 


i 


TR  8027 


»*» 

•S 


2w  G(c)  2  Boae)  -  B^e)  -  2g^  J^e)  - 


v, 

4’. 

.*,1 


—  Sr  Bq( re)  +  Qr  B]  ( re)  +  2gp  J.,(re)  + 


»*; 

»».. 

-ft 


+  4  -5  °n  B„  (ne)  - 

_  2  *  /  -»  D  0 


R  B.  (ne) 
n  1 


n=<P+2 


n=J?+2 


5 

l.l 


The  auxiliary  sequences  utilized  in  (30)  are  defined  below: 


SJL  gi+2  "  2g*+l  +  9A 


:1 

S 


S  =  g-2g.+g„ 

r  3r  r-1  r-2 


Qx  =JUX  *  1  )9j?+2  -  WJl  *  2)g^+]  +  (je  *  2)U  *  1  )9p 


3 


i*1 

5-1 


s*« 

£ 

:$ 

4, 


Qr  --  (r  -  2)(r  -  l)gr  -  2r(r  -  2)grl  4-  r(r  -  l)gp_2 


“n  ‘  V2  -  2Vl 


-  2Vl  -■  V? 


F  =  g  „  -  4g  ,-t-6g  -  4g  .  +  g  . 

n  yn+2  n+l  n  n-1  n-2 


R  =  n  0  +■  nF 

n  n  n 


for  n  - 

(i  *  2) ( 2) ( r  -  2) 


The  functions  B  (u)  and  B.(u)  are  those  defined  in  (13)-(17),  and 
o  1 

the  slash  on  the  summation  symbol  in  (30)  denotes  skipping  every  other  term, 


V’ 

■s 


after  starting  at  n  =  X  *■  2.  A  shorthand  notation  that  will  be  used  here  is 


n  -  X  +  2.  4  .  r  -  4,  r  -  2  =  (/  t-  2)(2)(r  -  2)  .  (33) 


r«j 

4*5 

»*5 

•f 
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Several  important  observations  should  be  made  about  the  result  in 
(30) -(32) .  The  four  quantities  in  (31)  are  evaluated  only  once  at  the  end 
points  n  =  JL  and  r.  The  sequences  in  (32)  must  be  evaluated  at  all  the 
points  listed  in  (33),  that  is,  at  every  other  interior  point.  All  of  these 
computations  should  be  done  once  and  stored,  when  given  the  function  g(x), 
the  limits  x^.x^,  and  sampling  increment  h,  prior  to  ever  considering 
which  cj  values  will  be  of  interest  in  (30).  Input  function  g(x)  must  be 
evaluated  at  all  x  =  xn  =  nh  for  n  =^(l)r. 

The  time-consuming  calculations  of  Bq ( u )  and  B^(u)  in  (30)  are  only 
necessary  at  the  values  u  =  ne  for  n  =y?(2)r,  and  need  not  be  evaluated  at 
any  of  the  in-between  points  n  =  ( +■  l)(2)(r  -  1).  The  Bessel  function 
J j  ( u )  need  only  be  evaluated  at  end  points  u  -■  and  re;  however,  this 
quantity  shows  up  as  a  free  by-product  of  evaluating  BQ(u)  and  B^u),  by 
the  method  indicated  in  appendix  A. 

SAMPLING  INCREMENT  FOR  u 

When  output  variable  <j  in  desired  integral  (27)  is  restricted  to 
multiples  of  a  sampling  increment  a,  according  to 

u  =  k A  for  1  £  K,  <  k  <  K2  ,  (34) 

then 

0  =  uh  -  kdh  ,  (35) 


and  (30)  takes  on  the  form 
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2kA  G(kA)  £ 


(kah) 


2  Sf  BQC£kAh)  -  B^kAh)  -  2g^  J-,  0?kAh) 


(kAh) 


-  Sf  BQ(rkAh)  Qr  B^rkAh)  +  2gr  J^rkAh)  + 


r-2 


r-2 


(kAh) 


H  °n  BQ(nkAh)  -  Rn  B^nkAh)  . 


(3b) 


n=J^-2 


n=J}+2 


At  this  point,  the  discussion  in  the  sequel  to  (26)  is  directly  relevant 
and  should  be  reviewed.  The  only  change  in  the  presentation  is  to  replace 


B^u),  there,  by  both  Bq(u)  and  B^(u)  here.  We  again  end  up  with  two 


alternatives  and  two  correspond ing  programs  for  evaluation  of  (36):  one 
faster  routine  which  may  require  considerable  storage,  and  a  slower  procedure 
utilizing  very  little  storage.  Programs  for  both  procedures  are  listed  in 
appendix  C. 


BEHAVIOR  FOR  SMALL  0 


When  e  is  small,  differences  of  functions  with  similar  values  are 


required  in  (30),  as  may  be  observed  by  the  linear  u>  dependence  on  the  left 
2 


side  and  the  1/0  dependence  on  the  right  side.  This  behavior  is  also 
typical  for  Filon's  method,  and  indicates  the  need  for  a  series  expansion  in 
powers  of  0  for  the  right-hand  side  of  (30)  when  @  is  small;  see 
[5;  ( 25 .4 . 53) J,  for  example.  The  appropriate  series  development  for  this 


parabolic  approx ima t i on  approach  to  (27)  is  given  in  (C-l  5) -(C-l  7) ,  through 


2  4  6 

order  0  .  Additional  terms  to  order  0  ,  0  can  be  derived  by  an 


obvious  extension  of  the  approach  given  there. 
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EXAMPLES 


Two  examples  will  be  considered  in  this  section;  the  first  is  a  Rayleigh 
function , 

g ( x )  =  x  exp(-x^/2)  for  x  >  0  ,  (37) 

for  which  8essel  transform  (11)  is  [9;  6.631  4] 

G(gi)  =  exp(-w2/2)  .  (38) 

The  second  is  a  Gaussian  function, 

2 

g(x)  =  exp(-x  )  for  x  >  0  ,  (39) 

leading  to  [9;  6.618  1  ] 

G(«)  =  l/2f?  exp(-u2/8)  Iq((j2/8)  .  (40) 

These  two  examples  are  very  different,  in  that  transform  (38)  decays  very 
quickly  for  large  u,  whereas  (40)  decays  very  slowly  for  large  a>.  In  fact, 
for  the  latter  case  [5;  9.7.1], 

G(w)  ~  l/«  as  w  -»  -t-o*  .  (41) 
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This  difference  will  enable  us  to  Investigate  both  absolute  and  relative 
errors  of  the  approximate  numerical  integration  procedures  developed 


earlier,  over  a  wide  range  of  values  of  u. 


ALIASING 


The  Bessel  function  JQ  is  rather  similar  to  a  sinusoid;  in  fact,  for 
large  z  [9;  9.2.1  ], 


•Uz)  ~  IT 


as  z  ■>  t*. 


(42) 


Then  when  argument  x  in  transform  (11)  is  sampled  at  increment  h,  we 
encounter  the  behavior 

Jo(a,Xn>  =  Jo("hn)  r  Jo(en)  ~  j  cos(en  -  4)  (43) 

for  large  en .  Now  when  e  -  2 ir,  the  cosine  yields  the  same  values  as  for 
0=0;  this  leads  us  to  expect  larger  errors  for  the  numerical  integration 
procedure  when  0  is  near  2ir. 


For  a  Fourier  transform,  this  aliasing  effect  was  studied  quant i tat i vely 
in  [10;  appendix  A]  for  both  the  Trapezoidal  rule  and  Simpson's  rule.  The 
former  rule  was  shown  to  have  a  large  aliasing  lobe  at  0  =  <*>h  =■  2«,  while 
the  latter  rule  had  an  additional  large  lobe  at  0  =  ir,  due  to  the 
alternating  character  of  the  Simpson  weights;  see  [10;  ( A  6 )  and  (A  -8 )  ] . 

This  leads  us  to  anticipate  that  the  linear  approx imat i on  procedure 
developed  here  for  Bessel  transform  (11)  will  be  subject  to  aliasing  near 
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8  =  2 *,  while  the  parabolic  approximation  will  be  degraded  earlier,  namely 
near  e  =  ir.  This  will  be  borne  out  by  the  numerical  examples  to  follow. 

GRAPHICAL  RESULTS 

The  Bessel  transform  numerical  integration  rule  for  the  linear 
approximation  to  g(x)  is  given  by  (23)  or  (26),  while  the  rule  for  the 
parabolic  approximation  to  g(x)  is  given  by  (30)  or  (36).  The  exact 
transforms  (38)  and  (40),  and  the  absolute  errors  associated  with  these  two 
rules,  are  depicted  in  figures  3  and  4  for  the  Rayleigh  and  Gaussian 
functions  g(x)  of  (37)  and  (39),  respectively,  with  sampling  increment 
h  =  .1.  The  ordinates  in  all  figures  are  the  logarithm  to  the  base  10  of 
the  correspond ing  results,  while  the  abscissas  are  linear  in  or  0.  The 
upper  limit,  xr,  of  integration  in  (2)  or  (11)  is  taken  large  enough  to 
guarantee  a  negligible  contribution  (less  than  IE-20)  to  the  truncation 
error. 

In  figure  3,  the  error  for  the  parabolic  fits  is  initially  lower  (for 
small  <j)  than  for  the  linear  fits;  however,  the  linear  error  decays  rapidly 
with  u),  and  stays  below  the  parabolic  error  for  larger  u.  Both  absolute 
errors  flatten  out  and  are  not  increasing  with  u,  at  least  for  this  range  of 
w  values.  The  maximum  value  of  0  is  .8,  as  indicated  in  the  figure. 


L0G10  Error  LOGlOlError 
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Figure  3.  Errors  for  Rayleigh  Function  g(x) 


Figure  4.  Errors  for  Gaussian  Function  g(x) 


For  the  Gaussian  function  g(x),  the  parabolic  error  in  figure  4  is 
everywhere  less  than  the  linear  error.  Both  errors  near  and  at  u  ^  0  are 
extremely  small;  this  fortuitous  result  for  the  linear  fits  is  fully 
explained  in  [6;  pages  92-93],  especially  in  the  paragraph  under  (3.4.5). 

It  has  to  do  with  the  fact  that  the  integrand  in  (11)  for  this  Gaussian 

2 

case,  namely  Jq(cjx)  exp( -x  ),  has  zero  odd  derivatives  at  the  limits  of 
integration.  This  is  not  the  case  for  the  Rayleigh  function;  hence  the  much 
larger  errors  at  u  =  0  in  figure  3  result. 

COMPARISON  OF  PROCEDURES 

To  demonstrate  the  benefits  to  be  accrued  from  the  fitting  procedures 
derived  in  this  study,  a  comparison  of  the  absolute  errors  for  four  different 
procedures  is  presented  in  figure  5  for  the  Rayleigh  function  (37).  The 
sampling  increment  in  x  is  h  =  .03.  The  variable  u>  now  covers  the  range 
(0,120);  the  point  where  8  =  ir  is  indicated  by  a  tic  mark  on  the  abscissa. 

The  Trapezoidal  result  is  obtained  by  applying  it  to  the  complete 
integrand  Jq(ujx)  g(x)  of  (11).  The  error  is  essentially  constant  for  all  u>, 
including  the  region  near  8  =  ir;  thus,  as  expected,  aliasing  is  not 
significant  at  8  =  ir  for  the  Trapezoidal  rule. 

Application  of  the  standard  Simpson's  rule  to  the  complete  integrand  of 
(11)  yields  a  very  small  error  near  w  -  0,  but  a  rapidly  increasing  error 
with  (j,  and  a  very  large  aliasing  lobe  centered  around  8  -  it.  Ihis  confirms 


the  expectations  presented  earlier  in  this  section. 
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Figure  5.  Errors  for  Four  Procedures;  £*-><120 


For  the  case  of  linear  fits  to  g(x),  rather  than  J0(ux)  g(x),  the 
error  drops  dramatically,  by  four  orders  of  magnitude  as  u  increases, 
similar  to  figure  3.  Furthermore,  there  is  no  aliasing  at  9  =  ir. 

The  situation  for  the  parabolic  fits  is  that  the  absolute  error  starts 
out  small  and  remains  so,  for  all  w  <  120,  there  being  a  slight  aliasing 
effect  near  8  =  ir.  However,  it  is  5  orders  of  magnitude  smaller  than  the 
Simpson  error  in  this  region  of  <e. 

The  results  in  figure  6  extend  the  abscissa  to  cover  the  range  of 
(120,240)  in  u;  that  is,  these  curves  are  an  extension  of  those  in  figure 
5.  Now  all  rules  suffer  aliasing  in  the  neighborhood  of  e  =  2ir.  The 
absolute  error  for  the  linear  procedure  increases  by  2  orders  of  magnitude 
near  e  =  2ir,  while  the  parabolic  error  is  just  slightly  larger;  however,  the 
latter  is  6  orders  of  magnitude  better  than  the  standard  Trapezoidal  and 
Simpson  rules  for  numerical  integration.  All  of  these  results  confirm  the 
predicted  presence  and  location  of  aliasing  discussed  earlier. 

ERROR  DEPENDENCE  ON  SAMPLING  INCREMENT 

In  figure  7,  we  investigate  the  dependence  of  the  error  on  increment  h 

employed  to  sample  x  in  (11).  Here  we  apply  the  linear  fit  procedure  to  the 

Rayleigh  function  (37).  The  absolute  error  for  small  u>  (<  2)  decreases  by  a 

factor  of  4  as  h  is  halved;  that  is,  the  large  error  bump  near  w  =  0 
2 

behaves  as  h  for  small  increments  h.  On  the  other  hand,  for  larger  ^  ( •>  5 ) , 


the  error  decreases  by  a  factor  of  16  when  h  is  halved;  that  is,  the 

4 

"saturation"  level  of  error  behaves  as  h  for  small  h.  The  slight  flare 
in  the  error  curve  near  <j  =  50,  for  h  =  .1,  is  an  indication  of  the 
beginning  of  aliasing;  that  is,  0  =  5  here,  which  is  near  the  0  =  2* 
location. 

Still  considering  the  Rayleigh  function  (37),  but  now  switching  to  the 
parabolic  procedure,  the  results  in  figure  8  demonstrate  that  the  error 

4 

drops  by  a  factor  of  16  as  h  is  halved;  thus,  the  error  dependence  is  h 
for  all  u.  The  wiggles  in  the  h  =  .1  curve  near  u  =  30  are  due  to  aliasing, 
since  0  =  tr  for  w  =  10*  =  31.4. 

When  the  function  g(x)  is  changed  to  the  Gaussian  example  of  (39),  and 

the  linear  fitting  procedure  is  employed,  the  errors  are  depicted  in  figure 

2 

9.  Here,  the  error  dependence  is  according  to  h  for  all  u>,  until 
aliasing  sets  in.  Aliasing  is  present  in  the  h  -  .1  curve  near  u  =  64, 
since  @  -  2*  at  u  =  62.8  for  that  curve.  Comparison  of  these  errors  with 
the  exact  answer  in  figure  4  reveals  that  the  relative  error  is  constant  in 
the  range  4  <  <j  <  56. 

When  the  parabolic  procedure  is  used  instead  on  the  Gaussian  example, 

4 

the  error  dependence  is  again  according  to  h  ,  until  aliasing  becomes 
dominant.  Ihe  aliasing  lobes  in  the  h  -  .1  curve  in  figure  10  are  centered 
at  0  =  *  and  2*,  as  before.  The  large  increase  in  the  error  for  the 
h  .025  curve,  when  u  exceeds  50,  is  a  feature  not  seen  previously.  It  may 
be  due  to  the  sum  of  distant  aliasing  of  sidelobes  which  decay  very  slowly 
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with  w;  in  fact,  from  (41),  the  exact  answer  only  decays  as  l/a>.  The  rapid 
decay  of  the  Rayleigh  transform,  (38),  apparently  precluded  this  type  of 
error  from  appearing  in  any  of  the  numerical  cases  considered  here  for  the 
Rayleigh  g(x) . 


28 
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SUMMARY 

T here  is  a  marked  difference  between  the  form  of  these  results  and  the 
Filon  equations;  namely,  the  term  multiplying  sample  value  gn  =  g(xR) 

(in  (B-3),  for  example)  varies  with  n  in  such  a  fashion  that  no 
simplification  or  factoring  is  possible.  In  order  to  better  explain  this 
complication,  let  us  investigate  the  evaluation  of  (18)  when  J^ux)  is 
replaced  by  exp(io>x);  that  is,  consider  evaluation  of  a  Fourier  transform, 
rather  than  a  Bessel  transform,  for  the  moment.  When  the  linear  fits  to 
g(x)  in  (18)  are  then  integrated,  there  follows 

ln  =  gn  h  exp(ine)  •  (44 

But  the  bracketed  term  here  is  a  common  factor  (independent  of  n)  that  can 
be  removed  from  the  summation  on  n.  This  fortuitous  simplification  does  not 
hold  for  the  corresponding  result  (B-3)  here,  because  whereas  exp(iu)  is 
periodic,  Jq(u)  and  A(u)  are  not. 

In  an  effort  to  recover  some  of  this  loss  in  execution  time,  we 
therefore  grouped  the  terms  in  (8-6)  in  an  alternative  form,  pivoted  around 
B^ne)  rather  than  gn;  see  (8-7).  Perhaps  another  rearrangement  of 
terms  would  be  more  advantageous  for  some  purposes. 
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It  is  possible  to  extend  the  results  here  to  other  Bessel  transforms. 
For  example,  suppose  we  are  interested  in  the  evaluation  of  first-order 
transform 


dx  J-|(o>x)  g(x)  , 


(45) 


and  we  approximate  g(x)  either  by  straight  lines  or  parabolas, 
integrals  in  (13)— (15)  are  then  replaced  by 

u 


The 


I 


dt  J1 (t)  =  l  -  jq(u)  , 


0 


u 

J  dt  t  J1( t)  =  B0(u)  , 


J  dt  t2  J1 ( t )  =  u2  J?(u)  -  2u  (u)  -  u2  Jq(u)  , 


(46) 


0 


where  we  used  [5;  (11.1.6)  and  (9.1.30)]  and  (16).  Since  all  of  these  terms 
have  already  been  encountered  here,  extension  to  transform  (45)  would  not  be 
difficult. 


For  the  evaluation  of  the  alternative  transform 


J ^  (  (jX  ) 

dx - - —  g( x )  , 


(4  7) 


30 
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we  need  the  additional  result  [5;  (11.1.1)] 


f 


Ji(t) 


dt 


4 

u 


I 


k  J2k(u)  = 
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=  u[J2(u)  +  2  J4(u)  +  3  J6^u)  4'--0  •  (48) 

But  this  type  of  term  is  easily  evaluated  by  means  of  the  downward 
recurrence  technique  given  in  appendix  A.  In  fact,  immediately  following 
the  single  line  Se  =  Se  +■  E,  we  have  merely  to  add  the  line  Sx  =  Sx  +•  Se; 
when  the  downward  recurrence  is  completed,  the  bracketed  term  in  (48) 
results  in  storage  location  Sx  (after  the  scaling  correction). 


31/32 
Reverse  Blank 
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APPENDIX  A 

NUMERICAL  EVALUATION  PROCEDURE  FOR  BESSEL  INTEGRALS 


The  three  fundamental  Bessel  integrals  that  must  be  evaluated  are  given 
by  (13) — ( 17)  as 

u 

A(u)  =  f  dt  JQ(t)  ,  (A-l) 

0 

u 

BQ(u)  =  A{ u)  -  u  JQ(u)  =  J  dt 

0 

u 

B 1  ( u )  =  A(u)  -  J1(u)  =  J  dt(l 

0 

By  expanding  in  a  power  series  [5;  (9.1.10)],  and  integrating  term  by 
term,  there  follows  from  (A-l), 


t  (u  -  t)  JQ(t)  ,  (A-2) 

"  u)  Jo(t)  ■  <#-3) 


A(u) 


a 


1 


2  ...k 

-u  / 4 ) 


=  u  - 


k=0  k !  k ! 


(k+j 


3 

u 

12 


u _ 

320 


( A-4 ) 


When  this  result  is  coupled  with  the  series  expansions  of  and  in 
(A-2)  and  (A-3)  respectively,  there  follows 


A-l 
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Bq(  u  ) 


t  and 

f 


) 

I 

I 


81  <u) 


(-u2/4  )k 

k!  (k  i-l)!  (k  4-  |) 


u 


80  2688 


u 

4 


V  (V/4  )* 

k=0  k!  ( k  -t-  1 )  !  (k  4-  i) 


3  5 

u  u_  _u _ 

2  "  48  *  1920 


( A  -  5 ) 


( A  -6) 


Although  these  power  series  could  be  used  for  small  and  moderate  values 

of  u,  they  are  not  useful  for  large  u,  due  to  the  loss  of  significant  digits 

caused  by  the  alternating  character  of  series  (A-4)-(A-6).  In  fact,  we 

will  find  that  a  downward  recurrence  will  yield  all  the  values  of  A,  B  , 

0 

B1 ,  J  ,  and  very  efficiently  for  small  u,  while  an  asymptotic  expansion 
is  equally  attractive  for  large  u. 


DOWNWARD  RECURRENCE 


We  start  with  [5;  (11.1.2)]  and  ( A - 1 ) : 

A ( u )  =  2[J,(u)  4  J  3 ( u )  4  J5(u)  4  . . . ]  .  ( A  -  7 ) 

Thus  if  we  can  evaluate  all  the  odd-order  Bessel  functions,  we  can  get  A(u) 

'  r')in  their  sum  Also,  8  (u)  and  8^(u)  follow  immediately  from  (A-2)  and 
>  ,  'I  we  i  an  idd'tiona’ly  get  J  (u). 


mi’  'fin  He  .  ,e  1  *  uni  T  arts  sat  >s*  y  t  he  lownwa  rd  ret  jr  -  eni  e  ;  S  ,  ■  4  .  '  .  ?  /  }  , 
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Jm<u>  W“>  -  Wu>  l*'8) 

for  m  >  0.  This  recurrence  can  be  started  by  guessing  at  JM(u)  =  0, 

J  (u)  =  IE-250  for  example,  and  evaluating  downward  via  (A-8)  to  id  -  0. 

M  ~  I 

Since  the  error  increases  much  slower  than  the  size  of  the  terms  in  (A-8) 

[5;  table  9.4],  the  relative  error  of  the  terms  is  very  small  for  the 
smaller  values  of  m,  if  M  is  chosen  large  enough  to  start  with.  In  order  to 
accurately  establish  the  absolute  level  of  the  sequence  of  [j  values,  we 
then  use  the  check  sum  formula  [5;  (9.1.46)] 

Jq(u)  +  2[J2(u)  +  J4(u)  +...]=  1  .  ( A —9 ) 

In  order  to  realize  15  decimal  accuracy  in  A,  B  ,  8,,  J  ,  J, ,  it  has  been 

o  1  o  1 

found  sufficient  to  choose  even  integer  M  as 

M  =  M(u)  =  2  I  NT  ^20  +  .  56u  -  jy  p'u)  +12  for  0  <  u  <  45  .  (A-10) 

While  conducting  the  downward  recurrence  on  m  in  (A-8),  an  even  sum  of 
Ju  f  i  +  •  •  •  .  and  an  odd  sum  of  Ju  .  <•  Ju  .  t  ...  ,  are  maintained. 

After  completion  to  m  =  0,  the  even  sum  is  subject  to  constraint  (A-9),  in 
order  to  establish  the  scale  factor  that  must  be  applied  to  all  the  desired 
outputs;  this  is  to  correct  for  the  initial  arbitrary  (incorrect)  guess  of 
Ju  ,(u)  -  IF. -250.  With  this  scale  factor  in  hand,  the  odd  sum  in  (A-7) 
can  then  be  modified  by  means  of  one  multiplication  for  the  correct  absolute 
level  for  A(u).  Since  the  last  two  quantities  yielded  by  recurrence  (A-8) 
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are  yu)  and  Jq( u )  (after  scaling),  we  then  have  all  the  necessary 
ingredients  to  determine  8q(u)  and  B^u). 

No  array  declarations  or  array  storage  is  necessary  in  this  procedure, 
since  there  is  never  any  need  to  "go  back  up"  the  recurrence  and  correctly 
scale  the  {yu)}  terms.  This  has  been  guaranteed  (through  numerical 
investigation)  by  the  choice  of  M  in  ( A— 1 0 ) .  A  further  economy  in  the 
program  for  this  two-term  recurrence  (A-8)  has  been  achieved  by  splitting  it 
into  even  and  odd  versions,  thereby  avoiding  the  usual  temporary  storage  of 
the  left-hand  side  of  (A-8)  until  the  right-hand  side  is  updated.  This 
compact  program  is  listed  below  as  subroutine  SUB  Bes j .  For  given  u,  it 
outputs  values  for  Jq(u),  J-j(u),  A(u),  Bq(u),  B^u),  provided  that  0  <  u  <  45. 

ASYMPTOTIC  EXPANSION 


For  large  u,  the  starting  integer  M  in  (A-10)  gets  too  large  to  make 
downward  recurrence  a  viable  procedure.  Instead,  we  resort  to  the 
asymptotic  expansion  [5;  (11.1.11)] 

u 

A(u)  =  |  dt  JQ(t)  ~  1  - 
0 


/  1 1  k 
(-1 )  a 

2k+l 

u 


2  k  +-1 


( A  —  1 1  ) 


as  u  ->  +  oo  ;  here,  we  also  used  the  definite  integral  result  that  A(co)  -  1 
[5;  (11.4.1/)].  The  values  of  the  coefficients  are  [5;  (11.1.2)] 


A -4 
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■  (*)„  i  (l)s 

s=0 

and  are  conveniently  obtained  by  recursion 

S"WS2SS!  25  Vl 

The  number  of  terms  required  in  the  summations  in  (A-ll)  depends  on  the 
value  of  u  and  the  desired  accuracy.  For  u  >  45  and  15  decimal  accuracy,  it 
has  been  sufficient  to  terminate  (A-ll)  at  k  =  INT(u/2). 

Since  (A-ll)  yields  only  A(u),  it  is  necessary  to  calculate  Jq( u )  and 
J  (u)  additional ly;  this  has  been  accomplished  by  use  of  [11;  section  6.8], 
All  of  these  quantities  are  evaluated  by  means  of  subroutine  SUB  Bessel 
listed  below.  For  input  u  >  0,  this  subroutine  yields  values  of  JQ(u), 

J] (u) ,  A(u),  Bo(u),  B] ( u) . 


2Ss! 


for  s  >  1 


(A-12) 


( A  -1  3) 


A  -5 


Ik  8U2 7 


10  SUE  Be  si e  1  '  X  ,  JO , J l , ft , B0, B 1  ■  1  ft  =  I  NTEGPftL  •  0  ,  X  ■  dr  Jo 

20  DOUBLE  k, I  '  INTEGERS 

30  IF  X  >45 .  THEN  60 

40  CALL  Besj <X, J0, J1 , ft, B0, B1  >  !  DGWNWftRD  RECURRENCE  9.1 


50 

SUBEXIT 

60 

I  =  I N  T <  X  )  /  2 

'  ftSYMPTOT I C 

SERIES  11.1 

70 

Rx=  1 .  /X 

30 

F  = , 5»Rx 

90 

T  = .  25 

100 

ft=1.25 

1  10 

Re  = . 625*Rx 

120 

I  m  =  P=  1 . 

130 

FOR  K= 1  TO  I 

140 

P  =  -P 

150 

S  ri  =  K  +  K 

1  60 

F  5  =  S  ri  -  .  5 

1  70 

F=F*F5*Rx 

1  SO 

T  =  T  *  F  5  x  k.  3  ri  +  S  n  > 

190 

ft  =  ft  +  T 

200 

Be  =  F*fl 

210 

Im= Im+P*Be 

220 

Sn=Sn+ 1 . 

230 

F5=Sn- . 5 

240 

F=F*F5*Rx 

250 

T  =  T*F5x(Sn  +  Sr"i> 

260 

ft  =  ft  +  T 

270 

Bo=F*fl 

280 

Re=Re+P*Bo 

290 

IF  Be*Be+Bo*Bo< 1 . E-26  THEN 

3  1  0 

300 

NEXT  K 

310 

F=X-. 78539816339744828 

320 

T=. 79788456080286541 

330 

ft  =  1  .  -  T  *  S  Q  R  <.  R  x  )  *  <  R  e  *  C  O  S  F  1  - 

I  rn  +  S  I  N  <  F  >  > 

340 

J0=FNJo<X) 

1  JO  =  J  O  \  X  > 

350 

J1=FNJ1(X> 

!  J  1  =  J  1  <  X  > 

360 

B  0  =  ft  -  X  *  J  0 

1  B  0  =  ft  (  X  > 

-  X.  J  0  '  X  > 

370 

B  1  =  ft  -  J  1 

'  B  1  =  ft  <  X  » 

-  J  1  '  X  > 

330 

SUEEND 

390 

i 

400 

SUB  Be i j  CU, J0, J 1 , ft, BO, El  > 

1  J  0  =  J  o  ■  U  1 

,  J  1  =  J  1  ■  u 

4  1  0 

IF  U>0.  THEN  450 

'  ft  =  ft  ■  U  ' 

=  I NTEGPftL  ■  O 

4  20 

J0  =  1  . 

i  B0  =  ft 1 U  < 

-  U  J  o  ■  U  > 

4  30 

Jl=ft=EO=Bl=0. 

1  B  1  =  ft  ■  U 

-  J  1  U 

440 

SUBEXIT 

450 

DOUBLE  Me, Ms 

!  INTEGERS 

.  27,  1 
11  &  12 


U  •  dt  J 


460  Mc=2*INT  <■  20.  + .  56*  U -175.  ■'  <  12.  +U  '  +  12 


470  T  =  2 .  U 


480 

Se  =  E  =  0 . 

490 

So=0= 1 . E -250 

500 

FOR  Ms=Mc  TO  2  STEP  -2 

510 

E  =  T*'s  Ms+1  >  * 0 - E 

'  9.1.27 

520 

S  e  =  S  e  +  E 

5  30 

0=T*Ms*E-0 

1  9.1.27 

540 

s  o  =  ■ j  o  +  0 

550 

NEXT  Ms 

560 

E  =  T  *  0  -  E 

570 

F  = l .  (Se+Se+E  > 

1  9.1.46 

580 

J0=E*F 

590 

J 1 =0*F 

6  0  O 

ft  =  >  S  o  +  S  o  >  *  F 

1  11.1.2 

6  1  0 

B0  =  ft-U*  JO 

620 

B  1  =ft-  J  1 

6  30 

SUBEND 

6  40 

i 

A-6 

650 

660 

670 

6  3  0 
6  '?•  0 

7  0  0 
7  1  0 
720 
730 
740 
750 
760 
770 
780 
790 
308 
810 
320 
830 
340 
350 
360 
370 
830 
390 
900 
9  1  0 
920 
9  30 
94  0 
950 
9  6  0 
-1  T-  0 
9  ;  0 
9  9  0 

1  0  0  0 
;  0 1 0 
1020 
1  0  3  0 
1040 
1050 
1  0  6  0 
10  7  0 
1  0  3  0 
1  0  9  0 
;  i  o  y 
1110 
1120 
1  1  3  0 

1  14  0 

1  1  5  0 
1160 
1  1  ~  0 
1  1  A' 


DEF  FNJo(X)  !  Jo(X)  via  Hart  #5845,  6546,  and  6946 

Y = R  B  3 ( X> 

IF  Y  >  8 .  THEN  770 
T  =  Y  *  Y 

F=227 1490439. 5536033-T+- 551 3534 . 564 770 7522- T -5292 . 6 1 7 1 30 33455 
P  =  2  3  34489 1 7 1 377869 . 7-T*<  4  776555944267  3 . 538-T+-  4621  72225031 . 71 
P= 1859623 1762 1397304. -T  + c 44 1 4532939 1 8 1 5932 . -T*P- 
0  =  204251433. 521  34357  +  T*  <.494030. 7949131  3972 +  T-*"  334. 720367561  75 
0=2344750013653996. 3  +  T  *  (  1581546244  9  769. 752  +  T**.  64398674535.  1  33 
0=185962317621897733. +  T  +  Q 
Jo=P/Q 
RETURN  Jo 
2  =  3.  '■< 

T  =  Z*Z 

Pn  =  2204. 5010439651304  +  T*<.1  23. 677535 7437141  9 +  T-.  908479 3 474 3 028 
Pn=8554 . 82254 1 50666 1 7  +  T* <  8894 . 4375329606 1 94  +  T *Pn  ' 

Pd  =  22 14. 04885191 47104  + T*< 138. 38490049992388+ T  ■ 

Pd  =  8554 . 32254 1 5866628  + T  +  *  3903. 836 1 4 1 7095954  +  T  +  Pd  > 

Qn=13. 990976865 960630  +  T*( 1 . 049 7327982 34554 8  + T* . 009352595  3294  0 
0  n  =  3  7 . 5105  349549571  12+T*<46. 0938268 1 4625 1 75  +  T*Qn > 

0 d  =  9 21. 5 6 697552653098 +  T  +  C74. 4 28 38974141  1  179  +  T  ■ 

Qd  =  2400 . 674237  1  1  72675  +  T*  c  297  1 . 9837452034920  +  T  *0d 
T=Y~ . 78539816339744323 

J  o  = .  282  09479177337820*  8  0  R  Z  '  *  C  0  3  *  T  >  *Pn  P  d  +  9  I  N  T  •  *  2  *  0  n  0  d  1 

RETURN  Jo 

FNEND 

BEF  F  N  J 1  <  X  >  1  J 1  ( X  >  Hart  #6045,  6747,  arid  7  147 

Y  =  R  B  3  >  X  ) 

IF  Y.3.  THEN  1040 
T  =  Y  *  Y 

P  =  .  1  1 07  35222445  37  306E- 1 0-T* .  6  3 1 94  3 1 0  3 1  744  3  1  6 1 E  -  1  4 
P= . 491 059  92  7655  5 1 294E-5-T ♦ •  . 9  382  1  9  3  365  1  407445E-3- T-P  • 

P  =  .  3  9  3  3  1  0  7  9  3  3  9  5  i!  3  3  -  T  *  '  .  1  7  0  5  7  6  9  2  6  4  3  4  9  6  1  7  1  E  -  2  -  T  »■  p  .> 

P  =  5373. 7877666  563200 -  T*',  61 . 2  t  3  7  6  9  9  7  3  5  6  9  4  3  9  -  T  *  P  1 
P  =  6  9  5  3  6  4  2  2 . 6  3  2  9  3  3  3  5  0  -  T  *  ■  8  3  5  6  7  8  5 . 4  8  7  3  4  3  9  1  4  3  -  T  *  ■  32  0  30  2.  74  6  3  8  5  3  9 
0=1  39072845. 265 9 676 9 +  T  +  * 6  7  0  5  3  4 . 6  3  3  5  4  8  2  2  9  9  3 ♦ T » ■  l  234 . 59  345  396 6 3 
J  i  =  *  p  -  u 

RETURN  J1 
2=3. ■  Y 
T  =  Z*Z 

P n  =  3  1  3  2.  7529 5 6  3 5 50695+T+,  1  74.  31  37974 3  379025  + T  *1.22  3  5 0 5  2:  7 6 4  3 5  9 
P  ri  =  1  2  9  0  9  .  1  3  4  7  1  3  9  6  1  3  3  1  +  T  *  i  1  3  0  '90. 4  2  0  5  1  1  0  3  5  0  6  5  ♦  T  »  F  n  • 

P d  =  3  1  0 9 . 2314  16  7  7 0 0 2 3 3 ♦  T *  ■.  169. 0 472177 5 0 0 3 6  1  0 ♦  T  • 

Pd= 1 2  909 .  1 3  4  7 1 396 1 879  +  T  *  <  l  3066 . 73 3087844020+ T +  Pd  1 
Or,  =  5  1  .  7  3 6 5  32 3  1  3  3 6 5 9  1  6  +  T  ♦  <  3  .  7994  4  5  3  7 9 6 9  3 0 6  7  3  +  T  ♦  .  0  2. 6  3 6  ;  4 6 6  4  760  2 

0  n  =  1  4  4 . 6  5  2  3  2  3  7  4  9  9  5  2  0  9  ♦  T  *  •  17  4, 4  2  9  1  6  3  9  0  9  2  4  2  5  9  +  T  »  0  r,  • 

0 d  = 1  1  1  9 .  1 0  9  3 5  2  7  0 4  7  4  3  7 ♦ T  * ■  35. 2  2  3 9  2 0 6  4  3  4  1  3  4  0 4  +  T 
'■  3=2  0 85 . 92  701  3  3  32  3  1  72  +  T  ♦  >  3  7  34.  34 0  1  0 6 0  1  6  301  •  ♦  T  * 0  j 
T*  i  -2.  3561944901  92  3448 

J  1  =  .  i  3  i  0  9  4  7  9  1  7  ^  i  8  7  8  2  0  •»  3  0  R  ■  Z  ■  ♦  1  L  Ci  '3  1  T  *  ♦  F  r,  F  .3  -  2  I  N  t  *  2  *  1  ’  r ,  J 

IF  0.  THEN  J 1 = - J  1 
REThRm  ji 
FuEMD 


74  ■ 

80  3-T-P  ■ 

504  +  T  >  1 

2  5  6  +  T  *  0  > 

30  3' 

3  1 9  > 


1 4  7  0  -  T  *  F 
;  0 1  9  +  T  ’  ' 

1 0  4  3  1 

■  4711  ■ 


A  -7/A-8 
Reverse  Blank 
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APPENDIX  B 

DERIVATION  OF  INTEGRATION  RULE  FOR  STRAIGHT  LINE  FITS  TO  g(x) 


The  situation  of  interest  here  is  represented  in  figure  1,  where 
straight  lines  are  fit  to  g(x)  between  adjacent  samples  of  g(x),  taken  at 
sample  points  hi-  In  particular,  the  contribution  to  integral  (11)  of 
an  (internal)  abutting  point  xn  was  set  up  in  (18)— (19) -  By  letting  t  =  ux 
in  (18),  and  using  (19),  (20),  and  (22),  namely 

x  -  x 

y  =  — — -  .  xn  =  nh  ,  e  =  o>h  ,  (B-l) 


there  follows,  for  the  n-th  contribution  to  the  integral, 

ne 

l  =  —  f  dt  J  (t)  (\  -  n  +  A  +• 
n  w  J  o  \  ©/ 

(n-l)e 
(nel  )e 

‘t;  I  dt  V')  (’  *n  -i)- 

ne 


( B-2 ) 


where  -  9(xn)  =  9(nh).  By  reference  to  the  auxiliary  functions  defined 
in  (13)-(17),  the  sum  of  these  two  integrals  can  be  expressed  in  the  compact 
form 
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The  procedures  in  appendix  A  are  now  directly  applicable  to  the  evaluation 
of  (B-3)  for  any  n. 

For  the  left-end  point  x^  depicted  on  the  left  side  of  figure  1,  the 
corresponding  contribution  to  desired  integral  result  (11)  is,  using  (B-l) 
again , 

Ijj  *  J  dx  JQ(wx)  g^  (1  -  y)  = 

(JM)e 

=  4-d  b1  Qje  4- 1)©]  -  1)  bJa©3  -  :),[>©]}  .  (8-4) 

The  corresponding  contribution  to  integral  (11)  for  the  right-end  point 
xr  is  given  by 

x 

r 

ir  =  §  dx  V-)  9r  f  y)  = 

x  -h 
r 

re 

■  Z  I  dt  V*> 

( r-1  )© 

g  * 

1  -  ])  BlT(r  -  U»]  -  (r  -  1)  B]  [re]  f  J,  [re]j  .  (B-5) 

B-2 


(As  a  check,  combination  of  (B-4)  and  (B-5),  upon  replacement  of  X  and  r  by 


n,  yields  (B-3),  as  it  should.  The  "end  correction  terms"  in  cancel 
out  for  all  internal  points,  n.) 


The  resultant  approximation  to  desired  integral  (11)  is  given  by  the  sum 


of  (B-3) -( B-5 ) : 


G(u) 


J  dx  J0(«*>x)  9(x)  35  Ig  +  Ip  +  ^  In  • 

X|  n=ft+l 


(B-6) 


This  particular  grouping  of  terms  is  according  to  the  function  sample  values 
=  [g(nh)}  .  An  alternative  grouping,  according,  to  the  samples  of 


function  B^u)  instead,  is  given  by 


’Pvi  ' 

a  < 

k  l)9x] 

B]  (Je) 

-  %  J1 

U0)  + 

*  V,  - 

(r  - 

•  l)9r] 

B]  ( re) 

*■  9r  J1 

(re)  + 

r-1 

♦  2. 

^n+1 

- 

f  Vi] 

B,(ne) 

n=|+l 

Whereas  B^ne)  must  be  evaluated  for  all  J  <  n  <  r,  the  function  need 
only  be  evaluated  at  the  end  points  ^0  and  re. 


TR  8027 


When  to  is  restricted  to  be  multiples  of  a  sampling  increment  A,  that  is 

id  =  kA  for  k  =  1  ,  2,  . . .  (B-8) 

then  (8-7)  yields,  for  k  >  1,  the  approximation 

ka  G(ka)  s  [P  g^+1  -  DgJ  B]  (JtkAh)  -  J^kAh)  * 

+•  [r  gr-1  -  (r  -  1  )gr[]  B1  ( rkAh)  -t-  g^  J^rkdh)  +• 

r-1 

+  2E_n[Vl  -  2gn  *  Vl]  Bl<nkah)  •  (B~9> 

n=^+l 

where 

gp  =  g(nh)  .  (B-10) 

Since  n  and  k  are  integers  (see  (20)  and  (B-8)),  the  evaluation  of 
B-|(u)  in  (B-9)  is  confined  to  integer  multiples  of  Ah,  i.e.  u  =  mAh. 

Further  discussion  on  how  to  take  advantage  of  this  feature  of  (B-9)  is 
given  in  the  sequel  to  (26).  The  end  result  is  that  we  have  two  alternative 
procedures  for  evaluation  of  (B-9)  and  two  corresponding  programs:  one 
faster  routine  which  may  require  considerable  storage,  and  a  slower 
procedure  utilizing  very  little  storage.  Programs  for  both  procedures  are 


listed  below. 


BEHAVIOR  FOR  SMALL  9 


When  e  is  small,  the  differences  of  like  quantities  in  (B-3)-(B-5)  can 
be  circumvented  by  expanding  B^  and  in  power  series  in  e.  Using  the 
facts  that 

3 

B-,(U)  ~  |  -  Jg  35  U  -»  0  , 

3 

Ji(u)  ~  f "  h  as  u  •* 0  *  (B-nj 


the  above  results  reduce  to 


3„h[l 

‘I  — 01- 


(8-12) 


as  9  -►  0.  By  use  of  the  power  series  expansion  developed  for  B^(u)  in 

4  6 

appendix  A,  these  results  could  be  extended  to  order  9,0  if  desired. 


The  total  contribution  to  (11)  is  given  by  the  sum  in  (B-6).  As  9  ■*  0, 
this  reduces  to  the  Trapezoidal  rule,  (12). 


10 

1  2ER0-TH  ORDER  BESSEL  TRANSFORM 

USING  LINEAR  INTERPOLATION 

20 

1  I NTEGRAL  <  X 1  ,  Xr  :>  dX  Jo(MX)  g<X) 

FOR  W 1 <  =W< =  l<12  IS  STORED  IN 

3  0 

1  G w •:  w h e re  N  =  K 5 * D €  1  w . 

F  as  t  e  r  h  i  g  h  -  s  t.  o  r  ag  e  . 

4  0 

De  l  .  =.025  i 

INCREMENT  rh  '  IN  X 

50 

L  =0  ! 

X 1  =  L  *  D  e  1  x  ,  L  •  a  0 

60 

R  =  4  0  0  ! 

X  r  =  R  *  B  e  1  x  ,  R  >  L 

?0 

D  e  1  w  =  .  2  1 

INCREMENT  (A;  IN  W 

80 

K  1  =  0  ! 

W 1 =K 1 *De 1 w ,  k  1  >  =  O 

80 

K  2  =  4  0  i 

W2  =  K2*De 1 w ,  K2  =K  1 

1  0  0 

DOUBLE  L , P , K 1 , K 2 . k.O , L 1 . R 1 , Hi, K 

s,I  '•  INTEGERS 

1  1  0 

D  I  M  G  1  5  0  O  1  ,  D  g  1  5  0  0  .>  ,  E  1  ■  5  O  U  O  O  ■  , 

J  1  1  '  1  0  O  >  ,  T  l  r  •  1  O  0  ■  ,  G  i.-i  1  1  O  0  ■ 

120 

10  =  1  1 

1  30 

K1=MRX< K  1,1* 

1  40 

L 1 =L  +  1 

150 

R 1 =R-  1 

1  6  O 

REDIM  Gx ( L : R • , Dg <' L 1 : R  > , B 1 < L*K  1 

:  R*K2> 

1  70 

R  E  D  I  M  J  1  1  <  k  1  :  k  2  <  ,  J  1  r  (  K  1 :  K  2  >  ,  G  w 

■;  k0 :  k2  > 

1  80 

FOR  K s  =  k 0  TO  K2 

180 

G  W  K  K.  S  >=0. 

200 

NEXT  ks 

210 

FOR  N i  =  L  TO  R 

220 

G> (Ns >  =  F  N  G  * N i *  D  e lx)  1 

SEE  DEF  FHG(X')  =  q  •  X  * 

2  30 

NEXT  Ns 

240 

G  1  =  G  >  ( L  > 

250 

Gr=G*.TR> 

260 

IF  K 0  > 0  THEN  320 

270 

F«.  5*  (Cl  *Gr  > 

280 

FOR  N s  =  L 1  TO  R 1 

2  80 

F  =  F  +  G  x  <  N  s  > 

3  0  0 

NEXT  Ns 

3  1  0 

G  w  i  0  j  -  F  *  D  e  1  . 

3  2  0 

FOR  N s  =  L  1  TO  R 

3  30 

D  g  c  M  s  >  =  G  <  x  N  s .  >  -  G  ■  N  s  -  1  ■ 

340 

NEXT  Ns 

350 

D2  =  Delw*D«  1 

360 

IF  L=0  THEN  410 

3  70 

FOR  k  =  =  k 1  TO  12 

3  8  0 

I =L*ks 

3  9  0 

CALL  Besse 1  1  I *D2,  JO,  J  1 1  ■  t  =  •  ,  A, 

B  O  ,  B  1  ■  I  ■  ■ 

400 

NEXT  Ks 

4  10 

FOR  ls=kl  TO  12 

420 

I  =  R  *  K  s 

4  3  0 

CALL  Bessel'  I  * D  2 ,  J O ,  J 1 r ■  1  s  •  , A , 

BO, B  1  '  I  ■  ■ 

440 

NEXT  Ks 

450 

FOR  N  s  =  L 1  TO  Rl 

4  60 

FOR  K  s  =  k  1  TO  1-2 

4  70 

I  =NS  i 

4  x  0 

IF  B 1  ■ -I  ■  0 .  Then  5O0 

4  **  i  j 

'.All  I  e  i  i  e  1  1  I  ♦  D  _  .  J  y  .  J  1  ,  A  .  E  u  ,  B  1 

,  I  .  , 

j  0  0 

NEXT  1  s 

/'J 


i '  *V * r  i- 
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Tl=L*Dg<Ll >-Gl 

T2=R*Dg  >'  R  >  -Gr 

IF  L  =  0  THEN  580 

F  OF  K  s  =  K  1  Tu  K  L 

G  w  1  K  s  '=TI»B1'L^  K  s  )  -  G  1  *  J  1  1  K  s  . 

NEXT  K s 

FOR  K  s  =  K  1  TO  K2 
F  =  T  2  *  B  1  <R*Ks>-Gr*Jlr  <  K.  s  ) 

G  ui  t  K  j)»Gui.Ks)-F 
NEXT  ks 

FOR  N  s  =  L  1  TO  R  1 
F  =  N  =  * 1  Bg1  Ns  +  1  •  -Bg1  Ns  ■ 

FOR'  k  s  =  k  1  TO  K  L 
G w  <  K s  ’>  —  G u  1  K  j  ,'  +  F*B  1  (.Ns* K s  ; 
NEXT  Ks 
NEXT  Ns 

FOR  K  s  =  K  1  TO  K2 
G w  <  K  s  >  =  G  w ( K  s ) / <  K s  * D  e 1 w> 

NEXT  Ks 
PRINT  Gw1*; 

PAUSE 

END 


D  E  F  F  N  G  i  X  > 

G  =X*EXPc-. 5*X*X; 
RETURN  Gx 
FHEND 


'  RAYLEIGH  EXAMPLE 
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16 
20 
30 
40 
50 
<£•  0 
70 
80 
90 
1  00 
1  1  0 
1  20 
1  30 
14U 
150 
1  8  0 
1  70 
1  80 

1  90 
200 
210 
220 

2  30 
240 
250 
280 
270 
280 
290 
300 

3  1  0 
3  2  0 
3  30 
-40 
350 
380 

3  70 
38  0 
390 

4  0  0 
4  1  0 
420 
4  30 
4  40 
4  5  0 
480 
4  70 
4  80 

4  9  0 
500 

5  1  0 


2ER0-TH  ORDER  BESSEL  TRANSFORM 
INTEGRAL (Xl  ,  Xr  >  dX  Jo<WX>  gOO 
Gw ( Ks  > ,  where  W  *  Ks*Delw. 

Del x  = . 025 
L  =  0 
R  =  4  0  0 
D  e  1  w  = .  5 
K  1  =0 
K  2  =  1  0  0 

DOUBLE  L , R , K 1 , K2 , K0 , L 1 , R  1  , N 
DIM  Gx<500) , Dg<500> , Gw<200> 
K0=K  1 

Kl=MAXua  ,  1  > 

L 1 =L  + 1 
R1=R-  1 


USING  LINEAR  INTERPOLATION 
FOR  W 1 < =W<  =W2  IS  STORED  IN 
Slower  low-storage. 
INCREMENT  <h>  IN  X 
X 1  =  L  *  D  e  1  x ,  L  >  =  0 
X  r  =  R  *  D  e  1  x  ,  R  >  L 
INCREMENT  IN  W 

Wl=Kl*Delw,  K 1 > = 0 
W2=K2*Del w,  K2>=K1 
Ks  1  INTEGERS 


RED  I M  Gx ’ L : R ) . Dg  <L1:R > , Gw  < K0: K2> 

FOR  K 3 =  K 0  TO  K2 
Lw(Ks)=0* 

NEXT  Ks 

FOR  Ns=L  TO  R 

G x  (.  N 3  j  =  F N G »'  N  s  *  D  e  1  x  )  !  SEE  DEF  FNGOO  = 

NEXT  Ns 
8  1  —  G  x  '■  l  ) 

Gr  =  GxC.  R) 

IF  K. 0  > 0  THEN  310 
F  =  .  5  *  t  G 1  +  G  r  > 

FOR  N  s  =  L  1  TO  R 1 
F  =  F  +  G  x  <  N  s  ) 

NEXT  Ns 

Gw<  0  >  =F*De 1 x 

FOR  N s  =  L 1  TO  R 

D  g  <■  N  *>«C  x  t  N  S  •  -  G  ...  (Ns-1  > 

NEXT  Ns 


D  2  =  D  e  1  w  *  D  e  1  x. 

T  l  =  L  ■*  B  g  <  LI > - G 1 
T  2  =  R  *  D  g  ■  R  >  -  G  r 
IF  L  =  0  THEN  430 
T  =  L  *  D  2 

FOR  K  s  =  K  1  TO  K  2 

CALL  Bessel-T*Ks,J0,Jl,A,B0,Bl  ■ 

GwU  S'=T1*B1-C1*J1 
NEXT  Ks 
T  =  R  *  D  2 

FOR  K  s  =  K  1  TO  K 2 

CALL  Bessel  ■  T  K  ;  ,  JO,  J  1  ,  A  ,  B  0  ,  B  1  • 
F  =  T2*B 1 -Gr * J 1 
G  w  1  s  *  =  G  w  v  K  s  >  -  F 
NEXT  Ks 

FOR  Ms  =  L 1  TO  R  1 
F  -Ns  *  ■  Dg  t  Ns  *  1  >  -  Dg  ■  Ns  ,■  ■ 

T=Ns-D2 


5  2  0  FOR  k s  =  1  1  TO  K  2 

5  10  CALL  Be  s  s  e 1 >  T -1 s  ,  T0 ,  j i , h , EO , B 1 

54  0  G  w  * K  s  ■  =  G  w 1 K  s  > * F  *  B 1 

550  NEXT  Ks 

580  NEXT  Ns 

5  70  F  0 R  (  s  =  f  l  JO  t  2 

5  :  M  ij'j'  >  s  1  =  Gw  I  ;  ■  I  s  *  Be  1  w  < 

5 90  NEXT  I  ; 


r  0>-'  PRINT  Gw* 

810  END 
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APPENDIX  C 


DERIVATION  OF  INTEGRATION  RULE  FOR  PARABOLIC  FITS  TO  g(x) 

The  situation  of  interest  here  is  depicted  in  figure  2,  where  parabolas 
are  fit  to  the  samples  of  g(x),  a  pair  of  adjacent  panels  at  a  time.  The 
derivation  of  the  resultant  approximation  to  integral  (27)  is  broken  down 
into  the  four  cases  illustrated  in  that  figure. 

It  is  again  presumed,  as  in  (20),  that. sample  points  of  g(x)  are  taken 
at  increment  h,  namely 

xn  =  nh  for  j(  <  n  <  r  ,  (C-l ) 

and  that,  in  addition, 

r  -J  is  even  .  ( C -2 ) 

That  is,  the  total  number  of  panels  employed  in  interval  x^  ,xr  must  be 
even.  A  breakdown  of  all  the  sample  points  |x^  into  the  four  categories 
of  figure  2  is  depicted  in  figure  C-l,  where  we  have  used  the  abbreviations 


LMAMAMAMAMR 


panel  panel 

pa i r  pa  i  r 


KM  r 


X 


Figure  C -1 . 


Ca t egor i /a t i on  of  Sample  Points 


v,y 

-'V 

vy 


,*,v, 

»*&* 

LI 

M 

•  ft 

m 

ir  ■ 


I 

m 

j 

s'V 

m 

$4 

i  jii 


M  =  mid-point  (of  a  panel  pair), 

A  =  abutting  point  (between  two  panel  pairs), 

L  =  left-end  point, 

R  =  right-end  point.  (C-3) 

It  is  presumed  in  the  following  that  a>  >  0;  the  case  for  o>  =  0  is  given  by 
(29),  while  u  <  0  is  immediately  covered  by  observing  that  Jq  is  even. 

Hid  -Point 

The  contribution  of  a  mid-point  xn  to  integral  (27)  is  (see  figure  2) 
x  +h 

V 

"n  '  J  dx  V“x>  9„  (1  -  V2)  » 

x  -h 
n 

(n+1 )e 

-9~Z  S  dt  Jo«l>  0.  "  (e  '  n)2]  •  (C"> 

( n-1 )e 

where  we  utilized  t  =  u>x,  ( C - 1 ) ,  (28),  and  (22).  Upon  expansion  of  the 


square  in  (C-4),  and  use  of  (13)-(17),  (C-4)  reduces  to  the  rather  compact 


( C — 5 ) 


Mn  =  ~2  f(n2  “  1}  t(n_1  )0]  -  B1  [  (n-Hl  )©] )  - 

-  -j  (Bo[(n-l)e]  -  Bo[(n+1)oj)}  . 

0 

This  type  of  term  is  yielded  for  n  =  j  t  1,  ^+3,  ...  ,  r-3,  r  -  1 ,  as 
reference  to  figure  C-l  will  verify.  Here,  and  in  the  following,  for  the 
sake  of  brevity,  we  do  not  document  the  rather  detailed  machinations  that 
lead  to  the  compact  form  (C-5)  from  the  integral  definition  (C-4).  The 
reader  will  have  to  reconstruct  those  nonprof itable  manipulations,  if 
interested. 

At  this  juncture,  instead  of  treating  an  abutting  point  with  its 
associated  4  panels  (see  upper  right  of  figure  2),  we  split  it  up  into  a 
panel  pair  with  a  left-end  point  and  another  panel  pair  with  a  right-end 
point.  We  thus  have  to  consider  a  general  left  point  and  a  general  right 
point. 

I.Ef-1  POINT 

This  case  is  obtained  by  looking  at  the  bottom-left  diagram  in  figure  2 
and  replacing,!  by  n  everywhere.  The  contribution  of  this  type  of  panel 
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v2h 

l.n  -  J  dx  J0(^x)  gn  (1  -  y)  (1  -  y/2)  - 
Xp 

(n+2)e 

I  «  J0(t)[!  -  f(|  -  n)  -  |(|  -  n)2]  - 
ne 

g 

=  2^{(n+1  )(n+2)(B1  [(n+-2)e]  -  B^ne])  - 

-  ^Bo[(n+2)e]  -  BQ[ne])  -  2  J^ne]}  .  (C-6) 

0 

This  type  of  term  is  yielded  for  n  =  i.  JU  2 . r  -  4,  r  -  2,  but  not 

n  =  r;  see  figure  C-l . 

RIGHT  POINT 


This  case  pertains  for  the  bottom-right  diagram  in  figure  2  when  r  is 
replaced  by  n  everywhere.  The  corresponding  contribution  to  integral  (27)  is 


r 


Rp  =  V  dx  JQ(ux)  gn  (1  f  y)  (1  +•  y/2)  = 


x  -2h 
n 


ne 


■  I  dt  V‘>  ['  *  l(i  - ")  -  |(i  - 


;n  -2)e 


y  „  r 

?^|(n-l)(n-2)(B1[ne]  -  B,  [  (n-2)e])  - 


4  (Bq[  ne]  -  Bo[(n-?)ej)  f  2  J,[ne]] 
0 


(C-7) 


C-4 


This  type  of  term  is  yielded  for  n  =  X  +  2,  Jl  +  4 . r  -  2,  r,  but  not 

n  -  X  ;  see  figure  C-l . 


ABUTTING  POINT 


We  can  now  immediately  obtain  the  integral  contribution  for  an  abutting 
point  (top-right  diagram  of  figure  2)  by  adding  (C-6)  and  ( C -7 ) : 

An  =  Ln  +  **n  = 

^n  f  1 

=  2^|(n+l )  ( n+2)  B1[(n^-2)0]  -  ~  BQ[(n^2)©]  -  6n  B^ne]  - 

0 

-  ( n-1 ) ( n-2 )  B1[(n-2)0]  +  BQ[ ( n-2)©]}  ,  (C-8) 

0 

which  holds  only  for  n  =  J  +■  2 ,  Jt  • t-  4,  ...  ,  r  -  4,  r  -  2;  see  figure  C-l. 

As  a  notational  shortcut,  we  say  n  -  (Jt  -t-  2 )  ( 2 )  ( r  -  2)  are  the  allowed 
values  of  n. 

At  this  stage,  we  have  succeeded  in  evaluating  all  the  types  of  terms 
that  have  been  depicted  in  figures  2  and  C-l.  The  total  approx imat i on  to 
integral  (27)  is  therefore 

r-1  r-2 

G(“>  s  S  \  <■  \  -  Rr  •  <C-9> 

n=jf+l  n=^+2 


in  terms  of  the  contributions  in  (C-5)-(C-8),  where  the  slash  on  the 
summation  symbol  denotes  skipping  every  other  term. 
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However,  this  grouping  of  terms  in  (C-9)  is  according  to  sample  values 
gn  =  g ( n h )  of  function  g ( x ) .  It  is  advantageous  to  re-arrange  this  sum, 
grouping  terms  instead  according  to  sample  values  of  functions  Bq(u)  and 
( u) ,  defined  in  (16)  and  (17).  After  considerable  manipulations,  the 
following  alternative  to  (C-9)  is  obtained: 

2<j  G(w)  a  B0C89)  -  B]  tfe)  -  2g^  J^e)  - 

0 

-  -j  Sr  Bq( re)  +  Qr  B] ( re)  +  2gp  J^re)  a- 

0 

r-2  r-2 

*4  V"S>  -  5-  R„  •  <C- 

0  n=j(+-2  n=4+2 

The  auxiliary  sequences  utilized  in  (C-10)  are  defined  below: 

\  =  gA*2  '  2Vl  "  9; 

Sr  =  gr  2Vl  +  gr-2 

Qf  *■  ng^2  -  na  *  2)qM]  f  (je  2>u  *  1^ 

0r  =  ( r  -  2 )  (  r  -  l)gr  -  2r(r  -  2)gr  }  *  r(r  -  1  )  9r^2  (c 

and 
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It  is  important  to  observe  from  (C-10)  that  the  Bessel  integrals  Bq(u) 
and  ( u )  need  be  evaluated  only  at  u  =  ne  for  n  =Jl(2)r,  and  need  not  be 
evaluated  at  the  in-between  points  n  =  +  l)(2)(r  -  1).  Of  course,  the 

input  function  g ( x )  must  be  evaluated  at  all  x  =  x^  =  nh  for  n  =^(l)r. 

The  quantities  in  (C-ll)  and  (C-12)  do  not  depend  on  9  ^  u>h,  and  can  be 
computed  just  once  and  stored,  in  preparation  for  use  in  (C-10). 

If  we  are  interested  in  evaluating  integral  G(o>)  in  (27)  at  values  of  u 
equal  to  integer  multiples  k  of  some  increment  A,  then  we  must  substitute 

o>  =  kA  and  e  =  u>h  =  kAh  ( C -1 3 ) 

into  (C-10).  Then  interest  centers  on  computation  of  BQ(u)  and  B^u)  at 
u  =  mAh  for  certain  integers  m.  This  consideration  has  been  discussed  in 
the  sequel  to  ( 26 ) . 

BEHAVIOR  FOR  SMALL  6 

When  e  is  small,  differences  of  functions  with  similar  values  are 
required  in  (C-10).  This  same  behavior  obtainsfor  Filon's  method;  see 
[5;  (25.4.53)].  Accordingly,  it  is  useful  to  have  a  series  expansion 
for  G((d)  about  e  =  0,  to  be  used  for  small  0 . 

Since  [5;  (9.1.12)] 

J0(u)  ~  1  -  \  u?  as  u  ->  0  ,  (C-14) 


C  -7 


•SSA6S&L 
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c 


substitution  in  (C-4),  along  with  the  change  of  variable  y  =  t/e  -  n,  yields 
the  mid-point  contribution 

1 

Mn  =  gn  h  y  dy  Jo(0(n+y))  (1  "  y2)  = 

-1 

~  gn  h  J  dy  [j  ■  4  ®2(n+y)2]  O  -  y2)  = 

-l 

=  3  9n  hf1  _  4  e?  (n?  +  5)}  as  9  0  •  ( c-l  5) 

A  similar  procedure  for  left  point  and  right  point  contributions  (C-6)  and 
(C-l)  gives 

L„  *  «„  ~  5  9„  h{l  -  4  eV  -  |>}  as  e  "*  0  .  (C-16) 

The  total  asymptotic  contribution  to  G(w)  in  (27)  is  therefore  given  by 
(a  modified  version  of  ( C -9 ) ) 

r-2  r  r-1 

g(cj)  ~  Ln  +  Rn  +  ^  Mp  as  e  ->  0  ,  (C-l/) 

n=Jt  n=J+2  n=4+l 

using  (C-l 5 )  and  (C-16).  For  9=0,  this  reduces  to  Simpson's  rule,  (29). 

4  6 

Additional  correction  terms  involving  9  ,  9  could  be  derived  by  using 
additional  terms  in  expansion  (C-14). 

When  gj  is  specialized  to  values  w  =  kA  in  ( C  - 1 0 )  ,  the  result  is  as  given 
in  (36).  Programs  for  both  a  faster  high-storage  procedure  and  a  slower 
low-storage  procedure  are  listed  below. 
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!  2ER0-TH  ORDEI 

20 

'  INTEGRAL  CXI 

3  0 

1  G  w  <  K  s  > ,  w  h  e  r ' 

4  0 

Be 1  <=. 03 

50 

L  =  0 

SO 

R  =  3  0  0 

70 

De  1  w=  1 . 

80 

K  1  =0 

90 

K  2  =  4  0 

1  00 

DOUBLE  L , R , K 

1  10 

DIM  G  x  i  8  O  0  ,  i 

120 

DIM  B 0 i 2  0000 

1  30 

K  0  =  K  1 

140 

K 1 =MAX  <  K  1  ,  1  > 

150 

L 1 =L+ 1 

160 

L2=L+2 

170 

R1=R-1 

180 

R2  =  R-2 

1  SO 

REDIM  GjULs  R 

200 

R  E  D  I  M  B  0  <  L  *  K 

210 

FOR  K s = K 0  TO 

220 

Gw<  Ks ) =0 . 

2  30 

NEXT  Ks 

240 

FOR  N  s  =  L  TO 

250 

Gx<Ns>=FNG<N: 

260 

NEXT  Ns 

270 

•  j  1  s  G  x  ^  L  } 

280 

G  r  -  G  ■  R  ■' 

2  9  O 

IF  K  0  >0  THEN 

3  0  0 

S  1  =  S  2  =  0  . 

310 

FOR  N s  =  L 1  TO 

Jlu 

S  1  =  S  1  G  x  >  N  s  > 

3  3o 

NEXT  Ns 

340 

FOR  N  s  =  L  2  TO 

3  5  0 

S  2  =  S  2  *  G  x  <.  N  s  > 

360 

NEXT  Ns 

3  70 

G  >o  1  O  >  -  ' G  1  +  G  r 

380 

1 J  1  1  =  1 J  A  1  L  l 

i  :*U 

G 1  2  =  G  - 1  L  2  > 

4  00 

G  r  1  =  G  <  (  R 1  > 

4  10 

ij  r  ^  =  U  a  v  R  2!  ■* 

4  20 

SI =G1 2-2. *G1 

4  30 

S  r  =  G  r  -  2 .  *  G  r  1 

440 

0 1 =L*L 1 *G 1 2- 

4  5  0 

Or =R2*R 1 *Gr - 

46  0 

G  1  2  =  G  1  *  2  . 

4  70 

G  r  2  =  G  r  *  2  . 

4  80 

D  d =  D  €  1  w  *  D  <c  1  ; 

4  90 

FOR  K  s  =  K  1  TO 

500 

F  =  K  s  *  D  2 

5  1  O 

b  1  K  s  '  =  t  •  '  1  F 

*5^0 

next  i; 

Xr)  dX  Jo<  WX>  g<X> 
e  W  =  K.  s  *  D  e  1  w . 


,  K  2 ,  K  0 ,  L  1  , 
i  w  <  5  0  @  >  ,  S  q  '■ 
,  B 1  '  20000  > 


USING  PARABOLIC  INTERPOLATION. 
FOR  W 1 < =  M <  =  W  2  IS  STORED  IN 
F  as  t  e  r  h  i  g  h  - s- 1.  o  r  ag  *  . 

INCREMENT  h  ■  IN  X 
X 1 =  L*De 1 x ,  L  =  0 
X r  =  R * D e- 1  x  ,  R-L  MUST  BE  EVEN  s« 
INCREMENT  (A)  IN  W 
Wl=Kl*Delw,  K 1 >  =  0 
1  W2=K2*De 1 w ,  K2>=K1 
L  2  ,  R  1  ,  R  2 ,  N  s  ,  K  s  ,  I  '  INTEGERS 
500  ■  ,  J  1  1  ■  5O0  ,  J  1  r  •:  500  1 


=  4 


:  K  1 : K  2  > ,  JlrCKl :K2> 


SEE  DEF  F N G  <  X >  =  g ( X 


;OU 


R 1  STEP  2 
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530  IF  L=0  THEN  580 


540 

FOR  k  s  =  k  1  TO  K2 

550 

I =  L  *  k  s 

5  8  0 

CHLL  Bessel  U*D2,  JO, 

Til'  Ks) , H , E 0 1  I  >,  B 

570 

NEXT  ki 

530 

FOP  k s  =  k  1  TO  k2 

530 

I  =  R  *  k  s 

8  0  0 

CHLL  Bessel  •:  I  *  B  2 ,  J  0 , 

J  1  r  <  k  S  '  ,  H  ,  B  0  I  >  ,  B 

810 

NEXT  ks 

820 

FOR  N s = L 2  TO  R2  STEP 

8  30 

FOR  k s  =  k 1  TO  k 2 

840 

I  =  N  s  *  k  s 

850 

IF  B0U9O0.  THEN  67 

0 

880 

CHLL  Bessel  U*B2,  J0, 

J 1 , H, B0(I > , B1 <  I  )  > 

870 

NEXT  ks 

880 

NEXT  Ns 

890 

IF  L=0  THEN  740 

700 

FOR  ks  =  Kl  TO  k2 

710 

I =L*ks 

7  20 

Gw < k  s ) »Sq < Ks >  *3 1  * B0  v 

I  .« -01  *Bl  <  I  >  -G1  2* J 

7  30 

NEXT  ks 

740 

FOR  ks=kl  TO  K2 

750 

I  =  R  *  k  s 

780 

F  =  8  q  (.  k  s  >  *  S  r  *  B  0  >:  I  :>  -Or 

*  B  1  <  I  >  -  G  f  2  *  J  1  r  i  k  s 

770 

Gw <  ks ) =0w <  Ks ) -F 

780 

NEXT  ks 

790 

FOR  Ns=L2  TO  R2  STEP 

800 

G  2  =  G  x  <  N  s  +  2  > 

3  1  0 

0  1  =  0  x  (.  N  S  +  1  } 

3  20 

H  1  =  G  x  <  N  s  -  1  > 

3  3  0 

H  eL  —  L  x  (  N  S  —  2  1 

34U 

B  n  =  G  2  -  2 . *  G 1 +  2 . *H1-H2 

8  50 

F  n  =  G  2  -  4  .  *  G  1  +  6  .  *  G  x  •-  N  s 

■-4. *H 1 +H2 

38  0 

R  n  =  N  s  *  <  N  s  *  B  n  ♦  F  n  > 

8  70 

FOR  k s  —  k 1  TO  k2 

380 

I  =  N  s  *  k  s 

3  9  0 

G  w 1  s  >  -  G  w  <  1  s  >  +  S  q  > 1  s  > 

-  E1  f -BO  I  - P' n -  E- 1  I 

900 

NEXT  ks 

1  0 

NEXT  Ns 

9  20 

F  =  B  e  1  w  *  2  . 

9  30 

FOP  ks=kl  TO  k2 

940 

IjuO  s^Gut  k S  '  1  1  s * F  1 

950 

NEXT  Ks 

980 

PRINT  G  w  <,  *  > 

9  70 

PH  USE 

980 

ENB 

9  9  0 

i 

1  0  0  0 

BEF  FNG'-X' 

'  g  ■  :  ■ 

1  0  1  0 

G  = .  >  E .  P  -  .  5  -  - 

1  PH-.LE 

1020 

RETURN  G  x 

1  O  3  0 

FNENB 

C-10 
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1  ZERO-TH  ORDER  BESSEL  TRANSFORM  USING  PARABOLIC 

INTERPOLATION 

20 

'  INTEGRALS  XI , X  r  >  dX  Jo<WX>  gOO  FOR  W1(=W<=W2  I 

3  STORED  IN 

30 

1  Gwf.Kj),  where  1*1  =  Ks*De1w.  Slower  low-sto 

a  '3  •=  • 

40 

De 1**. 03  i  INCREMENT  •  h * 

IN  X 

50 

L  =  0  !  X 1 =  L  *  D  e 1 x ,  L  >  = 

0 

r£.  0 

R  =  3  0  0  1  X  r =R*  D  e  1  :x  ,  R  -  L 

MUST  BE  EVEN 

?0 

De 1 w= 1  .  '  INCREMENT  (A; 

IN  U 

80 

K  1  =  0  !  W 1 =K 1 *De 1 w ,  K1 

>  =  0 

90 

K  2  =  1  2  0  1  U2=K2*Delw,  K2 

•  =  t  1 

1  0  0 

DOUBLE  L , R , K 1 , K  2 , K  0 , L 1 , L  2 , P 1 , R  2 , N s , K =  '  INTEGERS 

1  1  0 

D I M  G  1  800  1 , Gw ■  500 > , Sq 1  500  > 

1  20 

K  0  =  K  1 

1  30 

K 1 =MAX ( K  1  ,  1  ) 

1  40 

L 1 =L  +  1 

150 

L2=L+2 

1  80 

R 1 =R -  1 

1  70 

R2=R-2 

130 

RED  I  M  G r  L  :  R  > ,  Gw  <  K0:  K.2  ) ,  Sq  <  K  1 :  K2  3 

1  90 

FOR  K  s  =  K  0  TO  K  2 

200 

G  w  <  K  s  >  -  0 . 

210 

NEXT  K s 

220 

FOR  N  s  =  L  TO  R 

2  2  O 

G»‘Ni  1  =  FNG 1  Ni  *De  1  1  SEE  DEF  FNG’X1 

*  '3  '  '  ■  1 

2  40 

NEXT  Ns 

250 

U  1  =  G  x  1  L  .  * 

260  ■ 

G  r  =  G  ».  i.  R  1 

270 

IF  K 0 ; 0  THEN  360 

230 

S  1  =  S  2  =  0  . 

2  90 

FOR  N s  =  L 1  TO  PI  STEP  2 

O  0 

S  1  =  S  1  +  G  <  <.  Ns  > 

1 0 

NEXT  Ns 

3  _  0 

FOR  N s  =  L 2  TO  R2  STEP  2 

2;  2:  0 

S 2  =  S  2 ♦ G  a ■ N  s  * 

2  4  0 

NEXT  Ns 

2  5  0 

i  j  w  *  0  '  =  1 G  1  +  G  r  +  4  .  *  ;•  1  +  2  .  ♦  S  2  1  *  D  e  1  -  3  . 

2  6  0 

G  1  1  =  G  >•  >  L  1  ' 

7  0 

U  1  i  "  U  *  1  L  > 

2'  8  0 

G  r  1  =  G  *  R  1  > 

2  30 

Gr  2  =  G  •.  *  R 2  ■ 

4  00 

3  1  =  G  1  2  -  2  .  *  G  1  1  +  G  1 

4  1  0 

3  r  =  G  r  -  2  .  *  G  r  1  +  G  r  2 

4  2  0 

01=L*L1'*G12_2.'*'L^L2*G1  1+L2*L1*G1 

4  2  0 

Or  =  R  2  *  P  1  *  Gr  -  2  .  -*R* R2*Gr  1  +  R*R  1  »Gr  2 

440 

G  1  2  *  G  1  *  2 . 

4  50 

G  r  2  =  G  r  *  2  . 

4  6  0 

D  2  =  D  e 1 w*Bel 

4  70 

FOR  1  s  = 1  l  Tu  12 

4  20 

F  =  1  s  *  D  2 

4  9  0 

S  q  ■  Y  S'  =  1  .  •  F  »  F  * 

5  0  0 

NEXT  Is 

C-l  1 
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IF  L 
T  =  L  * 
FOR 
CALL 
G  w  <  K 
ME  XT 
T  =  R* 
FOR 
CALL 
F  = 

G  w 1 .  K 
NEXT 
FOR 
G2  =  G 
G  1  =G 
H  1  =G 
H2  =  G 
D  n  =  G 
F  r*i  —  G 
Rn  =  N 
T  =  Ns 
FOR 
CALL 
G  w K 
NEXT 
NEXT 
F  =  De 
FOR 
G  i.i 1  I 
NEXT 
FR  I  N 
FA  U 
END 


=3  THEN  570 

D2 

K. s  =  K  1  TO  K2 

Bessel  ■'  T  *  K  s ,  Jo,  J  1 ,  A ,  B  O ,  B  1  1 
s  >  =  S c| (  K s  '■>  * S  1  * B 0 - Q  1  *B1  - G  1  2 *  J  1 
Ks 
D2 

K  s  =  K  1  TO  K2 

Bessel  ■  T  *  K s ,  JO,  J 1 , A , BO , E  1  ■ 

1  K  s  >  *  S  r  ■*  B  0  -  0  r  *  E  1  -  G  r  2  *  J  1 
s  >  =  G  w  <  k  s  >  -  F 

•  K  s 

N s  =  L 2  TO  R2  STEP  2 
x  (  N  s  +  2  > 
x  C  N  s  ♦  1  > 
x  <  N  s  -  1  > 
x  <  Ns-2  > 

2-2. *G1 +2. *H1-H2 

2  -  4  .  *  G  1  +  6  .  *  G  x  N  s  ■  -  4  .  *  H  1  +  H  2 
s  *  i.  N  s  *  I>  n  +  F  n  > 

*D2 

Ks  =  K  1  TO  K.2 

Bessel ■  T *F  s  ,  JO ,  J 1 , A , BO , B  1  > 

S  >  *Gw <  ►'  s  >  +  Sq  1  K s  )  * Dn* B0-Rri* B  1 
K  s 
Ns 

1  w*2 . 

K S  =  K 1  TO  F  2 
S  1  =Gu  F  S  '  ‘  F  S  *  F  ) 

K  s 

T  Gui.  *  ' 


FNG*  > 

=  ;*E  :F"  5-x-x 

J  R  N  G  ^ 

ID 


'  9  ' 

1  F’  A  V  L  E  I  O  H  EXAMPLE 


•7:  -T-  jr  ’ 


’sT  TkT  V 


/-  l\  i-.  4  ^  .  - 


v  y  "$*  v 
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